%=====================
%this file solves the T2 model 
%=====================
clc;
clearvars;
global  climit  JR 
close all
%=====================
%define environment values
%=====================
JR=16; %retirement period 
r=0.10; %per period pre-tax interest rate
b=1; %no discount factor
%endowment for workers
e=[67
74
79
85
89
93
96
97
98
98
97
95
93
90
86
];

eJr=51; %endowment for retirees

climit=47.4; %utility consumption floor
s=9/10; %1/s is the expected lenghth for a constant survival rate s
amax=2500;
%====================
%Define state space and choice/value function
%=====================
NA=100; 
Avector=linspace(0,30,NA)';
Avector=Avector.^2; %more grids for small values
%for workers
vfunc=cell(JR,1); %the cell argument represents period
cfunc=cell(JR,1);
for j=1:JR
    vfunc{j}=zeros(NA,1);
    cfunc{j}=zeros(NA,1);
end %last period equals that of retirees
%for retirees
vrfunc=zeros(NA,1);
crfunc=zeros(NA,1);
%=====================
%Solve the model: recursively
%=====================

%start with retirees
critv=1.0e4;
diffv=1.0e-7;
crit=0.0004;
NCini=5; %initially check five points
while (critv>diffv)
    vrfunc_old=vrfunc;
          for ia=1:NA
            if ia==1
                cmin=climit;
            else
            cmin=crfunc(ia-1);
            end
            cmax=(eJr+(1+r)*Avector(ia)-tax_T2(eJr,r*Avector(ia)));
            cg=linspace(cmin,cmax,NCini);
            afuture=eJr+(1+r)*Avector(ia)-tax_T2(eJr,r*Avector(ia))-cg;
            vtemp=uc(cg)+b*s*interp1(Avector,vrfunc,min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=eJr+(1+r)*Avector(ia)-tax_T2(eJr,r*Avector(ia))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vrfunc,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=eJr+(1+r)*Avector(ia)-tax_T2(eJr,r*Avector(ia))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vrfunc,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            vrfunc(ia)=vmax_v;
            crfunc(ia)=cmax_v;
     
           end %ia
            
critv=max(abs(vrfunc_old-vrfunc));
end
%solve value function for workers
vfunc_temp=zeros(NA,1);
cfunc_temp=zeros(NA,1);
vfunc{JR}=vrfunc;
cfunc{JR}=crfunc;
for j=JR-1:-1:1 %retirees' life
    vfuture=vfunc{j+1};
        for ia=1:NA
            if ia==1
                cmin=climit;
            else
            cmin=cfunc_temp(ia-1);
            end
            cmax=(e(j)+(1+r)*Avector(ia)-tax_T2(e(j),r*Avector(ia)));
            cg=linspace(cmin,cmax,NCini);
            afuture=e(j)+(1+r)*Avector(ia)-tax_T2(e(j),r*Avector(ia))-cg;
            vtemp=uc(cg)+b*interp1(Avector,vfuture,min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=e(j)+(1+r)*Avector(ia)-tax_T2(e(j),r*Avector(ia))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=e(j)+(1+r)*Avector(ia)-tax_T2(e(j),r*Avector(ia))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture,min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            vfunc_temp(ia)=vmax_v;
            cfunc_temp(ia)=cmax_v;
            %vfunc{j}(ia)=vmax_v;
            %cfunc{j}(ia)=cmax_v;
     
        end %ia
        vfunc{j}=vfunc_temp;
        cfunc{j}=cfunc_temp;
end %j

%=====================
%simulate optimal decisions
%decisions are rounded to 1 dollar = 0.001 in thousands
%=====================
JN=100; %maximum periods
vopt=zeros(JN,1); %optimal value
copt=zeros(JN,1); %optimal consumption
aopt=zeros(JN+1,1); %optimal assets
taxopt=zeros(JN,1); %tax 
taxableopt=zeros(JN,1); %taxalbe income
%workers
for j=1:JR-1
            vfuture=vfunc{j+1};
            cmin=climit;
            cmax=(e(j)+(1+r)*aopt(j)-tax_T2(e(j),r*aopt(j)));
            cg=linspace(cmin,cmax,NCini);
            afuture=e(j)+(1+r)*aopt(j)-tax_T2(e(j),r*aopt(j))-cg;
            vtemp=uc(cg)+b*interp1(Avector,vfuture, min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=e(j)+(1+r)*aopt(j)-tax_T2(e(j),r*aopt(j))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=e(j)+(1+r)*aopt(j)-tax_T2(e(j),r*aopt(j))-cg;
                vtemp=uc(cg)+b*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            copt(j)=min(round(cmax_v,3),floor((e(j)+(1+r)*aopt(j)-tax_T2(e(j),r*aopt(j)))*1000)/1000);          

    aopt(j+1)=e(j)+(1+r)*aopt(j)-tax_T2(e(j),r*aopt(j))-copt(j);   
    taxopt(j)=tax_T2(e(j),r*aopt(j));
    taxableopt(j)=e(j)+r*aopt(j);
end
%retirees
for j=JR:JN 
            vfuture=vrfunc;
            cmin=climit;
            cmax=(eJr+(1+r)*aopt(j)-tax_T2(eJr,r*aopt(j)));
             %start with three points
            cg=linspace(cmin,cmax,NCini);
            afuture=eJr+(1+r)*aopt(j)-tax_T2(eJr,r*aopt(j))-cg;
            vtemp=uc(cg)+b*s*interp1(Avector,vfuture, min(amax,afuture),'spline');
            [vtemp_v, ctemp_g]=max(vtemp);
            cmax_v=cg(ctemp_g);
            vmax_v=vtemp_v;
            
            diff=(cmax-cmin)/4;
            while diff>crit

            diff=diff/2;
          
            if ((cmax_v>cmin))
                cg=cmax_v-diff;
                afuture=eJr+(1+r)*aopt(j)-tax_T2(eJr,r*aopt(j))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                    continue
                end %if
            end %while internal
            if ((cmax_v<cmax))
                cg=cmax_v+diff;
                afuture=eJr+(1+r)*aopt(j)-tax_T2(eJr,r*aopt(j))-cg;
                vtemp=uc(cg)+b*s*interp1(Avector,vfuture, min(amax,afuture),'spline');
                if (vtemp>vmax_v) 
                    vmax_v=vtemp;
                    cmax_v=cg;
                end %if
            end %while internal
            end %while
            copt(j)=min(round(cmax_v,3),floor((eJr+(1+r)*aopt(j)-tax_T2(eJr,r*aopt(j)))*1000)/1000);          
    aopt(j+1)=eJr+(1+r)*aopt(j)-tax_T2(eJr,r*aopt(j))-copt(j);   
    taxopt(j)=tax_T2(eJr,r*aopt(j));
    taxableopt(j)=eJr+r*aopt(j);
end


%payment from the experiment if adopting the optimal strategy
paymentopt=zeros(JN,1);
j=1;
paymentopt(j)=b^(j-1)*uc(copt(j));
for j=2:JN
paymentopt(j)=paymentopt(j-1)+b^(j-1)*uc(copt(j));
end




%payment if adopting the hand to mouth strategy
paymenthandtomouth=zeros(JN,1);
j=1;
paymenthandtomouth(j)=b^(j-1)*uc(e(j)-tax_T2(e(j),0));
for j=2:JR-1
paymenthandtomouth(j)=paymenthandtomouth(j-1)+b^(j-1)*uc(e(j)-tax_T2(e(j),0));
end
for j=JR:JN
paymenthandtomouth(j)=paymenthandtomouth(j-1)+b^(j-1)*uc(eJr-tax_T2(eJr,0));
end

chandtomouth=zeros(JN,1);
j=1;
chandtomouth(j)=e(j)-tax_T2(e(j),0);
for j=2:JR-1
chandtomouth(j)=(e(j)-tax_T2(e(j),0));
end
for j=JR:JN
chandtomouth(j)=(eJr-tax_T2(eJr,0));
end

%the probability of decease in period j
prob_d=zeros(JN+1,1); 
for j=JR+1:JN
    prob_d(j)=prob_d(j-1)+(1-s)*(1-prob_d(j-1));
end
prob_d(JN+1)=1; %maximum age of JN
fprintf('Expected tax revenue')
sum((1-prob_d(1:JN)).*taxopt)



%rename optimal decision with a T2 tag
copt_T2=copt;
paymentopt_T2=paymentopt;
taxableopt_T2=taxableopt;
taxopt_T2=taxopt;
sopt_T2=zeros(JN,1);
aopt_T2=aopt;


for j=1:JR-1
sopt_T2(j)=e(j)-taxopt_T2(j)-copt_T2(j);
end

for j=JR:JN
sopt_T2(j)=eJr-taxopt_T2(j)-copt_T2(j);    
end


fprintf('T2: Selected aggregate statistics')
[sum((prob_d(2:JN+1)-prob_d(1:JN)).*paymentopt_T2); paymentopt_T2(15); sum((prob_d(2:JN+1)-prob_d(1:JN)).*paymentopt_T2)-paymentopt_T2(15);...
    0; aopt_T2(JR);0;  mean(copt_T2(1:JR-1));  mean(copt_T2(16:25)); copt_T2(16); sum((1-prob_d(1:JN)).*taxopt_T2); sum(sopt_T2(1:JR-1))]

%save the key variables
save T2opt.mat copt_T2 paymentopt_T2 sopt_T2  aopt_T2 taxableopt_T2 taxopt_T2   prob_d


fprintf('T2: decision rules\n')
fprintf('1.consumption   2.cumulative earnings   3.saving    4.asset balance  5.ordinary income   6.taxes  7.survival prob.')
T2=[copt_T2(1:JN) paymentopt_T2(1:JN) sopt_T2(1:JN)  aopt_T2(1:JN) taxableopt_T2(1:JN) taxopt_T2(1:JN) (1-prob_d(1:JN))]







